Exact Current Statistics of the ASEP with Open Boundaries 
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Non-equilibrium systems are often characterized by the transport of some quantity at a macro- 
scopic scale, such as, for instance, a current of particles through a wire. The Asymmetric Simple 
Exclusion Process (ASEP) is a paradigm for non-equilibrium transport that is amenable to exact 
analytical solution. In the present work, we determine the full statistics of the current in the finite 
size open ASEP for all values of the parameters. Our exact analytical results are checked against 
numerical calculations using DMRG techniques. 
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A system containing carriers (of thermal energy, mass, 
or electrical charge) and subject to a driving field in its 
bulk, or to a difference of potentials between its bound- 
aries, will usually evolve to a non-equilibrium steady 
state (NESS) with a non-vanishing macroscopic current 
of heat, particles or charges flowing through it. Due to 
the presence of this macroscopic current, time-reversal 
invariance is violated. This is a situation which lies be- 
yond the realm of traditional thermodynamics: steady- 
state transitions at the microscopic level break detailed 
balance and the principles of equilibrium statistical me- 
chanics do not apply. Hence, for a system that is bulk- 
driven, boundary-driven, or both, no suitable generaliza- 
tion of the Gibbs-Boltzmann formalism exists that would 
allow us to predict the value of the current and of its fluc- 
tuations from first principles. 

During the last two decades, substantial progress has 
been made towards a statistical theory of non-equilibrium 
systems [THE]. Large deviation functions, that encode 
atypical fluctuations of a physical observable, are likely 
to be the best candidates to generalize the traditional 
thermodynamic potentials. Moreover, it has been proved 
that large deviations functions display symmetry prop- 
erties, called 'Fluctuation Theorems', that remain valid 
far from equilibrium [2]. These remarkable relations im- 
ply linear response theory in the vicinity of equilibrium. 
Hence, the determination of large deviations in a non- 
equilibrium system, whether theoretically, numerically, 
or experimentally, is a question of fundamental impor- 
tance [MI]. 

There are very few models in non-equilibrium physics 
that can be studied analytically. Among these, the asym- 
metric simple exclusion process (ASEP) has become a 
paradigm [T5lll7j . The ASEP is a one-dimensional lat- 
tice gas model in which particles perform biased random 
walks and interact through an exclusion constraint that 
mimics a hard-core repulsion: two particles cannot oc- 
cupy the same site at a given time. This minimal sys- 
tem appears as a building block in a great variety of 
phenomena that involve low-dimensional transport with 
constraints. Invented originally to represent the motion 
of ribosomes along mRNA during protein synthesis, this 
model plays a seminal role in non-equilibrium statistical 



mechanics and has been applied to problems as different 
as surface growth, biological transport, traffic flow and 
pure mathematics [5j [T8H22] . 

In the long time limit, the ASEP reaches a NESS with 
a fluctuating macroscopic current. Exact results have 
been derived for the exclusion process on a periodic ring 
and on the infinite line, using the Bethe Ansatz, deter- 
minantal processes and random matrix theory |22H26j . 
For open boundaries, the steady-state has a recursive 
structure [57] that can be encoded by a matrix prod- 
uct representation [2S] , a fruitful method to analyze low- 
dimensional transport models [29] . The mean value of 
the stationary current, the associated density profiles and 
the phase diagram of the open ASEP are known exactly 
[2jj[28]. However, finding the full statistical properties of 
the current in the open ASEP has remained, until now, 
an outstanding challenge that has stimulated many works 
[51 [TT1 15DH55] . A recent conjecture based on the Bethe 
Ansatz [30] gives the asymptotic behavior of the large 
deviation function of the current for infinitely large sys- 
tems in some specific regions of the phase diagram. In the 
present work, we give exact analytic expressions for the 
full statistics of the current, that are valid for arbitrary 
system sizes and boundary parameters, thus solving this 
long-standing problem. 

The dynamics of the ASEP is that of a continuous 
time Markov chain: during an infinitesimal time interval 
dt, a particle located on a site can jump forward to the 
next adjacent site with rate 1 and hop backward to the 
previous site with rate q, provided these sites are empty. 
A particle can enter the site 1 with rate a and the site 
L with rate <5, and can exit from the site 1 with rate 7 
and from the site L with rate (3 (see Fig. [lj. Each of 
the 2 L microscopic configurations C of the ASEP can be 
written as a binary string of length L, (ti, . . . , t£), where 
Tj = 1 if the site i is occupied and Tj — otherwise. The 
probability Pt(C) of being in configuration C at time t 
evolves according to the Master equation: 

dP ^=Y J M{C,C')P t {C). (1) 
a 

The non-diagonal matrix element M(C,C) represents 
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the transition rate from C to C. The diagonal element 
M(C,C) = - YjC=^c M(C',C) is equal to minus the exit 
rate from C. 
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FIG. 1. Dynamical rules for the ASEP with open boundaries. 
The rate of forward jumps has been normalized to 1. Back- 
ward jumps occur with rate q < 1. All other parameters are 
arbitrary. 



In the long time limit, the ASEP reaches a NESS where 
each configuration C occurs with a probability P*(C), 
that can be written as a matrix product |28j : 

1 L 

P*(C) = r Wll (nD + (l - n)E) \v) , (2) 

2—1 

where the operators D and E, the bra-vector (W\ and 
the ket-vector \V) satisfy quadratic algebraic relations 



EE- qED = (1 - q) (D + 
((3D-5E)\V) = (l-q)\V) 
(W\(aE~jD) = (l-q)(W\ . 



E) 



(3) 



The normalization constant in equation ^ is given by 
Z L = (W\ (D + E) L \V) . The matrix product represen- 
tation allows us to determine stationary equal-time cor- 
relations and density profiles for any system size L. 

For L — > oo, the ASEP has three phases whose bound- 
aries are given in terms of the effective densities p a — 
l/(a+ + 1) and pb = b + /(b + + 1) of the left and right 
reservoirs, where 



(1 - q - a + 7) ± - q - a + -f) 2 + 4cry 
fl± = Ya ' (4) 



, (1 - g - 13 + S) ± y/(l - g - I + 8) 2 + 4/M 

b± = 2? ' (5) 

The ASEP is in the Maximal Current (MC) phase when 
p a > 1/2 and pb < 1/2, in the Low Density (LD) phase 
when p a < 1/2 and p a + Pb < 1 and in the High Density 
(HD) phase when pb > 1/2 and p a + Pb > 1 • 

For a system of size L, the average value J of the 
stationary current is given by the ratio Z^-i/Z^ which 
can be expressed in terms of orthogonal polynomials |34j . 
However, the fluctuations of the steady-state current can 
not be calculated from the knowledge of the stationary 
probabilities alone. In order to study the current we in- 
troduce an observable Y t that counts the number of parti- 
cles exchanged between the system and the left reservoir 



between times and t. Therefore, Y t +dt — Y t + y where 
y = +1 if a particle enters the site 1, y = —1 if a parti- 
cle exits from 1 during the interval dt and y = other- 
wise. These three mutually exclusive types of transitions 
lead to a three parts decomposition of the generator M: 
M = M + + M- + Mq. We note that Y t also represents 
the total integrated current that has flown through the 
system till time t. When t — > 00, the expectation value 
of Y t /t converges to the average stationary current J. 
The convergence rate is quantified by the large deviation 
function <I>(j), characterizing non- typical fluctuations of 
Y t and defined as P = j) ^e' m J\ 

A different manner to encode the statistics of Y t is 
through its characteristic function which, in the long 
time limit, behaves as (e^ Yt ) ~ e £ ^) 1 ; where £(p) is 
the cumulant generating function of Y t , and is the Leg- 
endre transform of the large-deviation function [14) : 
£(p) = maxj (pj — <&(j)). Following OE], one can prove 
that £(p) is the largest eigenvalue of the deformed opera- 
tor M(p) = e^M + +e~"M_ +M . Thus, the calculation 
of the cumulants of the current is equivalent to an eigen- 
value problem. 

For the ASEP with periodic boundary conditions, 
M{p) can be diagonalized by Bethe Ansatz, leading to 
a full solution for the current fluctuations O El]. In 
the case of open boundary conditions, intcgrability con- 
ditions are only met on hypersurfaces of the parameter 
space and the Bethe Ansatz can be used only for L — > 00 
and in specific regions of the phase diagram |30 . 

We have obtained a solution valid for all parameters 
and all system sizes using a generalized matrix product 
representation. The components F^(C) of the dominant 
eigenvector F a of M (p) can be expanded formally as a 
power-series with respect to p to any given order k > 0. 
For each value of k, we have proved rigorously [35] that 
F^ can be represented by a matrix product Ansatz up to 
corrections of order p k+1 i.e., 



1 L 

F ^°) = —) ( W k\ II ^ Dk + \ V k)+° 
%L i=l 

. (6) 

The matrices E^ and Ek are constructed recursively 
starting with E\—E and E\= E and 

D k+1 = (1 ® 1 + d ® e) ® E k + (1 ® d + d ® 1) ® E k 

E k+X = (1 ® 1 + e ® d) <8> E k + (e ® 1 + 1 ® e) ® D k 

(7) 

where we have defined the operators d = D — 1 and e = 
E — 1 that satisfy the q-deformed harmonic oscillator 
algebra de — qed = 1 — q. These matrices are related 
to the ones used for the matrix product solution of the 
multispecies periodic ASEP [36] . 

The boundary vectors (W k \ and |Vfc) are constructed 
by taking tensor products of bra and ket vectors. We 
start with \Vi) = \V) and (Wi\ = (W\ and iterate 

|14+i) = \V) ® \V) ® \V k ) (8) 
(W k+1 \ = {W»\®(W'>\®(W k \, (9) 
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where \V) is defined in Eq. ^ and 

[0(l-d)-*(l-e)]|V)=O 



(10) 



(W|[a(l + e"e) -7(l + e^d)] = (l-g)(W|(ll) 



(W"|[a(l-e? e)-7(l 



; d)]=0. 



(12) 



This matrix Ansatz allows us to calculate the cumu- 
lants to any desired order k. Our central result is a 
parametric formula for the cumulant generating function 
£(»); 



ST, 



k>l 



and £ = — (1 — q) D k 



(13) 



fc>i 



where -B is a formal parameter that has to be eliminated 
from the two equations. We emphasize that similar para- 
metric expressions have appeared in all known exact ex- 
pressions for the current cumulant generating function 
0[23J|37] and a similar generic form was derived from the 
additivity principle in [S]. The function £(p) is fully spec- 
ified from the knowledge of the scalars Ck and Dk- These 
are given by contours integrals in the complex plane along 
a contour T (to be defined below): 



C k 



dz 4>k(z) 
2iir z 



Di 



dz 4>k(z) 
r 2i^(z + iy 



(14) 



The 4>k(z)'s are obtained as follows: we define a function 
Wb(z) that depends on the parameter B 



k>l 



(15) 



and we find that Wb(z) is uniquely determined as the 
solution of the functional equation: 



W B {z) 



ilnll 



BF{z)e 



X[W B ]{z) 



(16) 



where F{z) is given by the expression 



(l+*) L (l + *-T(* J )oo(*-'%o 

(a+z)oo(^ t )oo(a-^)oo(^)oo {b + z) 00 ( h ^) QO {b_z) 00 { h f) 00 

(17) 

with (x)oo = rifclo(l ~ 1 kx ) ■ We n0 ^ e * na * F( z ) appears 
in the definition of the Askey- Wilson polynomials, known 
to be relevant to the open ASEP [33] ■ The operator X 
is a linear integral operator: 



X[W B ](zi) = / 7 ^W B {z 2 )K 
J T 2^7^ z 2 

where the kernel K is given by 

00 k 



z k + z k 



} 



(18) 



(19) 



k=l 



and the contour T in the complex plane encircles (once) 
the points 0, q k a + ,q k a^,q k b + and q k b- for all integers 



k > 0. The kernel K(z\j z 2 ) was used in [23] to calculate 
the current fluctuations in the periodic case. 



The functional equation ( 16 ) contains the complete in 



formation about the current statistics: by solving it iter- 
atively to any order k, we obtain the first k cumulants of 
the current. At first order, we have </>i(z) = F(z)/2 and 
the mean value of the current is 



J=hm^ 

t— ¥00 t 



(1-9) 



Di 



(1-?)- 



£_ 

rr 2i7 



F(z) 



dz F{z) 
JT 2ztt (z+1) 2 



(20) 



This expression is identical to that given in [34] . At sec- 
ond order, the variance of the current is: 



where C 2 and D 2 are obtained using (14) with 



i 



F 2 (z) + 



dz 2 F{z)F{z 2 )K{z/z 2 ) 

2lTTZ 2 



For higher cumulants, exact expressions similar to 
Eq. (21) are obtained and can be expressed via a com- 



binatorial tree expansion akin to that found in the peri- 
odic case [24]. The expression of the diffusion constant 
A generalizes the formula of [3S] obtained for the to- 
tally asymmetric exclusion process (TASEP) in which 
q = 6 = 7 = 0. For the TASEP, the kernel K and 
the operator X vanish identically and F(z) reduces to 



Ftasep(z) 



(1 + z) 2L (l- z 2 ) 



2\2 



z L (l - az)(z - o)(l - bz){z - b) 



(22) 



and b — 



_ 1-/3. 



then, Eq. (16) leads to 



with a = - 

4>k{z) = -Ftasep(- z )/2 an d the results of [37J are retrieved. 
For a periodic system of size L with N particles, the 
current fluctuations can be brought into the framework 
described here with the same generalized matrix Ansatz, 
but the boundary vectors are replaced by a trace and 
equation (16) is modified as follows: the prefactor 1/2 
is removed, F(z) = (1 + z) L / z N and the Kernel K is 
still given by (19). Then, the results of [23j, originally 



obtained by Bethe Ansatz, are retrieved. 

The derivation of the above results involves combina- 
torial identities for matrix elements of the generalized 
matrix Ansatz. Some of these identities were guessed by 
induction rather than mathematically proved [38] . It was 
therefore necessary to validate our calculations numeri- 
cally. For small size systems (L < 7), expressions for the 
cumulants have been checked against the exact values 
from direct calculations. For larger systems (L < 100), 
we compared the analytical formulas with numerical com- 
putations of the cumulants performed using a DMRG 
method. That method, originally introduced to study 
ground state properties of quantum spin chains |39j was 
recently adapted to calculate the highest eigenvalues of 
deformed stochastic operators like M(fi) [TTJ. A few of 
those results are displayed in Figs [2] and [3] 
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FIG. 2. Third (E3, red) and fourth (E4, blue) cumulants in 
the Maximal Current phase, with q = 0.5, a+ = b+ = 0.65, 
a_ = b_ = 0.6; the full lines represent the corresponding large 
size asymptotic behaviors. 
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FIG. 3. Second (E2, red) and third (E3, blue) cumulants in 
the High Density phase, with q = 0.5, a+ = 0.28, 6+ = 1.15, 
a- = -0.48 and = -0.27. 



For large system sizes, L — > oo, the cumulant gen- 
erating function takes different expressions in different 
phases. These are derived from an asymptotic analysis 
involving the leading singularities of F(z) [3H [33] • In 
the Low Density phase, the dominant singularity is the 
pole at a + leading to <j>k{z) ~ F k {z). Using the Lagrange 
inversion formula as in [37], we obtain 



£(A0 = (l-g)(l-p„) 



e** + (l- Pa)/ Pa 



(23) 



This expression agrees with the one found in [3U] using 
the Bethe Ansatz. Its Legendre Transform matches the 



prediction of the additivity principle 
$0') = (1 - (?) \ Pa - r + r(l - r) In 



^ Pa 
Pa 1 



(24) 

-r). 

> b. 



where the current j is parametrized as j = (1 — q)r(l 
The HD Phase leads to similar expressions with a + - 
and p a — > 1 — p . In both cases, the statistics of the 
current do not depend on system size in the large L limit. 

In the Maximal Current phase, we find that fc-th cu- 
mulant grows as 

i (*-3)/2_ when L 

— > oo, we have 

-1/2 



£ - 



p = 



E 



(2fe)! 



^(k+3/2) 



(1 ~ q )L 
16 



l 

-3/2 



fe=l 



(2fc)! 



(25) 



B k .(26) 



These expressions have the same structure as those ob- 
tained for the case of a periodic ring [9 and the large 
deviation functions have the same asymptotic behavior. 
Moreover, for the open TASEP of size L with a = 1 and 
ft = 1/2, we observed that the formulas are identical to 
those for the half-filled periodic TASEP of size 2L + 2. 
Along the shock line (p a = 1 — < 1/2), we obtain 



p = -2L 



-i (! + «+) 



OO 



2fc-l 



(l-a H 



(2fc)! 



(27) 



£ 



(1 - q )a + 
(l + o 



P 



-2LT 



: (1 ~ g)a- 



E 

fe=i 



,2fe-2 



(2fc) 



(28) 



with the fc-th cumulant scaling as L^'~ 2 ^ as can be ex- 
plained by the domain wall picture for p a « 1 [32l [37] . 
We note that this is the only case where the statistics of 
the current depend on both the system size and on the 
boundary parameters at the large L limit. 

We have obtained exact formulas for the current statis- 
tics of the open exclusion process in contact with two 
reservoirs. Our results are valid for arbitrary sizes and 
values of the parameters and have been tested by precise 
DMRG computations in various regions of the phase di- 
agram. They could also be used as benchmarks to test 
alternative computational algorithms [12 . In the limit of 
large size systems, the asymptotic behavior of the large 
deviation function is derived in all regions of the phase 
diagram as long as the asymmetry (1 — q) remains fi- 
nite. The diffusive limit q — > 1 represents an important 
open analytical problem and the exact formulas should 
coincide with the predictions of macroscopic fluctuation 
theory [7J [31] • We have used an extension of the ma- 
trix Ansatz, that was introduced for multispecies exclu- 
sion models [36j . The relation between multispecies mod- 
els and current fluctuations (and also between open and 
periodic systems) is mysterious as no direct mapping is 
known. We believe that our results should be derivable 
from Bethe Ansatz for a spin chain with non-diagonal 
boundaries, but the corresponding Bethe equations have 
not yet been derived [30 . Finally, we emphasize that the 
matrix representation given here contains all the infor- 
mation about the density profiles that generate atypical 
currents: the precise calculation of these profiles repre- 
sents a challenging open question [T]. 
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